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ABSTRACT 

We investigate the mass profile of ACDM halos using a suite of numerical simulations 
spanning five decades in halo mass, from dwarf galaxies to rich galaxy clusters. These 
halos typically have a few million particles within the virial radius (7200), allowing 
robust mass profile estimates down to radii below 1% of r2oo- Our analysis confirms 
the proposal of Navarro, Frenk & White (NFW) that the shape of ACDM halo mass 
profiles differs strongly from a power law and depends little on mass. The logarithmic 
slope of the spherically-averaged density profile, as measured by f3 — — dlnp/dlnr, 
decreases monotonically towards the center and becomes shallower than isothermal 
(j3 < 2) inside a characteristic radius, r_2- The fitting formula proposed by NFW 
provides a reasonably good approximation to the density and circular velocity profiles 
of individual halos; circular velocities typically deviate from best NFW fits by less 
than 10% over the radial range which is well resolved numerically. On the other hand, 
systematic deviations from the best NFW fits are also noticeable. Inside r_2, the 
profile of simulated halos becomes shallower with radius more gradually than predicted 
and, as a result, NFW fits tend to underestimate the dark matter density in these 
regions. This discrepancy has been interpreted as indicating a steeply divergent cusp 
with asymptotic inner slope, 0q = (3{r — 0) ~ 1.5. Our results suggest a different 
interpretation. We use the density and enclosed mass at our innermost resolved radii 
to place strong constraints on (3q: density cusps as steep as r -15 are inconsistent with 
most of our simulations, although /3q = 1 is still consistent with our data. Our density 
profiles show no sign of converging to a well-defined asymptotic inner power law. We 
propose a simple formula that reproduces the radial dependence of the slope better 
than the NFW profile, and so may minimize errors when extrapolating our results 
inward to radii not yet reliably probed by numerical simulations. 



1 INTRODUCTION 



Disk galaxy rotation curves; strong gravitational lensing 
by galaxies and clusters; the dynamics of stars in ellipti- 
cal galaxies and of gas and galaxies in clusters; these are 
just examples of the various luminous tracers that probe 
the inner structure of dark matter halos. Such observations 
place strong constraints on the distribution of dark matter 
in these highly non-linear regions that may be contrasted di- 
rectly with theoretical predictions for halo structure. Given 
the sensitivity of such predictions to the nature of the dark 
matter, these observational constraints constitute provide 
tests that may question or even rule out particular models 
of dark matter. 

Robust prediction of the inner structure of cold dark 
matter (CDM) halos is a quintessential N-body problem, al- 



beit one of considerable complexity due to the large overden- 
sities and, consequently, the short crossing times involved. 
Indeed, only recently have computational capabilities im- 
proved to the point of allowing realistic simulation of the 
regions which house the luminous components of individual 
galaxies. 



This work builds upon the pioneering efforts of Frenk 
et al (1988), Dubinski and Carlberg (1991), and Crone et 
al (1994), among others, which led to the identification of 
a number of key features of the structure of dark matter 
halos assembled by hierarchical clustering. One important 
result of this early work concerns the absence of a well de- 
fined central "core" of constant density in virialized CDM 
halos. In this sense, dark matter halos are "cuspy": the dark 
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matter density increases apparently without bounds toward 
the center of the halo. 

A second result concerns the remarkable similarity 
("universality") in the structure of dark matter halos of 
widely different mass. This was first proposed by Navarro, 
Frenk & White (1996, 1997; hereafter NFW), who suggested 
a simple fitting formula to describe the spherically-averaged 
density profiles of dark matter halos, 

p ^ = i i \rTT7 i m > W 
(r/r s ){l + {r/r 3 )Y 

where p 3 and r a are a characteristic density and radius, re- 
spectively. The larger the mass of a halo, the lower its char- 
acteristic density, reflecting the lower density of the universe 
at the (later) assembly time of more massive systems. 

Further simulation work of similar numerical resolution 
(see, e.g., Cole & Lacey 1996, Huss et al 1999) provided 
support for the NFW conclusions, although small but sys- 
tematic differences began to emerge as the numerical resolu- 
tion of the simulations improved (see, e.g., Moore et al 1999, 
hereafter M99, Ghigna et al 2000, Fukushige & Makino 1997, 
2001, 2003). These authors reported deviations from eq. 
that increase systematically inward, and thus are particu- 
larly noticeable in high-resolution simulations. In particular, 
Fukushige & Makino (2001) reported that NFW fits tend to 
underestimate the dark matter density within the scale ra- 
dius r s . M99 reached a similar conclusion and interpreted 
this result as indicating a density cusp steeper than that of 
the NFW profile. These authors preferred a modified fitting 
function which diverges as r" 1 ' 5 near the center, 



PW (r/r M ) 1 - 5 ( 1 + (.r/rM) 1 - 5 )' w 

One should note, however, that there is no consensus 
amongst N-body practitioners for such modified profile (see, 
e.g., Klypin et al 2001 and Power et al 2003, hereafter P03), 
and that there is some work in the literature suggesting that 
the central density cusp might actually be shallower than 
r" 1 (Subramanian et al 2000; Taylor & Navarro 2001; Ri- 
cotti 2003). 

This unsettled state of affairs illustrates the difficulties 
associated with simulating the innermost structure of CDM 
halos in a reliable and reproducible manner. The high den- 
sity of dark matter in such regions demands large numbers 
of particles and fine time resolution, pushing to the limit 
even the largest supercomputers available at present. As a 
result, many of the studies mentioned above are either of in- 
adequate resolution to be conclusive or are based on results 
from a handful of simulations where computational cost pre- 
cludes a systematic assessment of numerical convergence. 

Establishing the detailed properties of the central den- 
sity cusp, as well as deriving the value of its asymptotic 
central slope, are important for a number of reasons. For 
example, steeper cusps place larger amounts of dark matter 
at the center, exacerbating the disagreement with observa- 
tions that suggest the presence of a constant density core in 
low surface brightness galaxies or in strongly barred galaxies 
(Moore 1994, Flores & Primack 1994, McGaugh & de Blok 
1998, Debattista and Sellwood 1998, van den Bosch et al 
2000). Steep cusps would also be important for direct detec- 
tion experiments for dark matter, as a possible gamma-ray 
annihilation signal of WIMPS at the Galactic center would 



be particularly strong for r cusps (Calcaneo-Roldan & 
Moore 2000; Taylor & Silk 2003; Stoehr et al 2002). 

Finally, the detailed structure of the central cusp is not 
the only focus of contrasting claims in the literature. For ex- 
ample, the "universality" of CDM halo structure has been 
questioned by Jing & Suto (2000), who find that the slope 
of the density profile at a fixed fraction of the virial radius 
steepens towards lower halo masses. Klypin et al (2001), on 
the other hand, point out that such a systematic trend is 
entirely consistent with universality as originally claimed by 
NFW, and just reflects the mass dependence of halo char- 
acteristic density. 

We address these conflicting issues here using a suite 
of nineteen high-resolution simulations of the formation of 
halos in the standard ACDM cosmogony. Halo masses are 
chosen in three main groups: "dwarf" halos with M200 ~ 
10 10 h,- 1 M Q , "galaxy" halos with Af 2 oo ~ 10 12 h' 1 Mq and 
"cluster" halos of mass M200 ~ 1O 15 /i _1 M . This allows 
us to gain insight into the effects of cosmic variance at each 
mass scale, as well as to explore the mass dependence of the 
structure of ACDM halos. We define the mass of a halo to 
be that contained within its virial radius, that is, within a 
sphere of mean density contrast 200 1 . 

This paper is organized as follows. Section [5] describes 
briefly the numerical simulations; § |3] discusses our main 
results; and we summarize our conclusions in §2] 



2 NUMERICAL EXPERIMENTS 

The numerical set up of our simulations follows closely the 
procedure described by P03, where the interested reader 
may find a thorough discussion of our initial conditions gen- 
crating scheme, the choice of N-body codes and integrators, 
as well as the criteria adopted to optimize the choice of the 
numerical parameters of the simulations. For completeness, 
we include here a brief discussion of the main numerical 
issues, but refer the reader to P03 for a more detailed dis- 
cussion. 

2.1 N-body codes 

The simulations reported in this paper have been performed 
using two parallel N-body codes: GADGET, written by Volker 
Springel (Springel, Yoshida & White 2001), and PKDGRAV, 
written by Joachim Stadel and Thomas Quinn (Stadel 2001). 
As discussed in P03, both codes give approximately the same 
results for appropriate choices of numerical parameters, and 
neither code seems obviously to outperform the other when 
similar numerical convergence criteria are met. 

2.2 Cosmological Model 

We adopt a flat, fio = 0.3 ACDM cosmological model whose 
dynamics is dominated at present by a cosmological con- 
stant, £7a = 0.7. The matter power spectrum is normalized 

1 We use the term 'density contrast' to denote densities expressed 
in units of the critical density for closure, p cl -it = 3H 2 /8nG. We 
express the present value of Hubble's constant as H(z = 0) = 
H = 100 h km s" 1 Mpc" 1 
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so that the present linear rms amplitude of mass fluctua- 
tions in spheres of radius 8 /i -1 Mpc is erg = 0.9. We assume 
a linear fluctuation power spectrum given by the product of 
the square of the appropriate CDM transfer function, T 2 (k), 
and a Harrison-Zel'dovich primordial power spectrum (i.e. 
P(k) ock). 

2.3 Parent Simulations 

The halo samples were drawn from three different ACDM 
cosmological "parent" simulations. Table Q lists the main 
numerical parameters of each of these simulations: Lbox is 
the size of the cosmological box, Zi is the initial redshift, 
m p is the particle mass, and e is the softening parameter, 
assumed fixed in comoving coordinates. 

The dwarf, cluster, and most of the galaxy halos, were 
extracted from the simulations ACDM-512 (Yoshida, Sheth 
& Diaferio 2001) and SGIF-128. These two parent simula- 
tions, both carried out within the Virgo Consortium, used 
the CDM transfer function given by CMBFAST (Seljak & 
Zaldarriaga 1996), assuming h — 0.7 and Qt = 0.04. This 
transfer function is well fit by the Bardeen et al (1986) 
fitting formula with a value of 0.17 for the shape param- 
eter r. Three of the galaxy halos (G1-G3, see Table |2J 
were extracted from the parent simulation described by Eke, 
Navarro & Steinmetz (2001, labeled ENS01 in TableHJ. That 
simulation used the Bardeen et al (1986) fitting formula for 
the CDM transfer function, with h = 0.65 and T = 0.2. 

2.4 Initial Conditions 

Since completing the numerical convergence tests reported 
in P03, we have developed a more flexible and powerful set of 
codes for setting up the resimulation initial conditions. This 
resimulation software enables us to iterate the procedure in 
order to "resimulate a resimulation", an important step for 
setting up appropriate initial conditions for dwarf halos. The 
basic methodology employed is very similar to the methods 
described in P03, with just a few minor differences. Galaxy 
halos G1-G3 were selected from the ENS01 simulation and 
their initial conditions were created using the software de- 
scribed in P03. All of the other halos were set up with the 
new codes, following the procedure we describe below. 

The first stage is to carry out, up to the redshift of in- 
terest (typically z = 0), a "parent" simulation of a large, 
representative volume of a ACDM universe. These parent 
simulations are used to select halos targeted for resimula- 
tion at higher resolution. Once a halo has been selected for 
resimulation at z = 0, we trace all particles within a sphere 
of radius ~ 3 V2oo to the z = oo "unperturbed" configura- 
tion. We then create a set of initial conditions with much 
higher mass resolution in the volume occupied by the halo 
particles, and resample the remainder of the periodic box at 
lower resolution, taking care to retain sufficient resolution 
in the regions surrounding the halo of interest so that ex- 
ternal tidal forces acting on the high-resolution region are 
adequately represented. 

The procedure involves two main steps. Firstly, we set 
up a uniform multi-mass distribution of particles to approx- 
imate the particle positions in the high-resolution region at 
z — oo. This is accomplished by arranging particles either in 



a cubic grid or as a "glass" , within a cube just big enough to 
contain the region of interest. Either choice approximates a 
uniform mass distribution very accurately. Outside the cube 
we lay down particles on a set of concentric cubic shells, cen- 
tered on the cube, which extend outwards until they fill the 
entire periodic volume of the parent simulation. These con- 
centric shells are filled with more massive particles whose in- 
terparticle separation increases approximately linearly with 
distance from the high resolution region. Unlike the grid or 
glass, this arrangement does not reproduce a perfectly uni- 
form mass distribution. However, by populating each shell 
with regularly spaced particles, we obtain a configuration 
which is uniform enough for our purposes. 

In the interest of efficiency, we replace those particles 
in the high-resolution cube that do not end up in the se- 
lected halo with more massive particles made by combining 
several high-resolution ones. This procedure, particularly for 
the dwarf halo resimulations, significantly reduces the num- 
ber of particles in the initial conditions and the run time of 
the subsequent simulation. Thus, each halo forms from an 
"amoeba-shaped" region consisting only of the highest reso- 
lution particles in the hierarchy. We have explicitly checked 
that the resampling procedure adds no extra power; in tests, 
the multi-mass particle distribution remains very close to 
uniform over an expansion factor of up to ~ 50. 

Once a multi-mass but uniform mass distribution has 
been created, the next step is to to add the appropriate 
Gaussian density fluctuations. This is done by assigning a 
displacement and a peculiar velocity to each particle using 
Fourier methods. By using the same amplitude and phase for 
every Fourier mode present in the parent simulation, a per- 
turbed density field essentially identical to that of the parent 
simulation can be reproduced. In the high-resolution cube, 
because the particle mass is smaller than in the parent sim- 
ulation, it is necessary to add additional short wavelength 
modes (with amplitudes fixed by the adopted power spec- 
trum) down to the Nyquist wavelength of the new particle 
grid. To ensure that the Fourier transforms needed to add 
this extra power are of a manageable size, we make the addi- 
tional power periodic on the scale of the central cube rather 
than on the scale of the parent simulation. The longest wave- 
length added is typically smaller than one tenth of the side 
length of the original cube. As described in P03, the indi- 
vidual components of the displacement field are generated 
in turn, and the displacements calculated at the particle po- 
sitions by trilinear interpolation. To set up growing modes, 
we use the Zel'dovich approximation and make the peculiar 
velocities proportional to the displacements. 

The initial redshift, Zi, of each resimulation is chosen 
so that density fluctuations in the high-resolution region 
are in the linear regime. P03 find that convergent results 
are obtained when Zi is high enough that the (theoretical) 
rms mass fluctuation on the smallest resolved mass scale, 
a(m p , Zi) does not exceed ~ 0.3 (where m p is the mass of a 
high- resolution particle). All of our simulations satisfy this 
criterion. 

2.5 Halo selection 

The resimulated halos analyzed in this paper were all iden- 
tified in the parent simulations by applying the friends-of- 
friends (FoF) group finding algorithm (Davis et al 1985) 
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with a linking length I = 0.164. Cluster-sized halos were 
drawn from a 479 ft -1 Mpc simulation volume (ACDM-512 
in Table 0. The FoF(0.164) groups were first ordered by 
mass and then ten consecutive entries on the list centered 
around a mass of 10 ft -1 Mq were selected. Galaxy halos 
were likewise drawn from a 35.325 /i -1 Mpc volume (SGIF- 
128), with the exception of three of the halos (G1-G3) which 
were selected from a 32.5 h -1 Mpc volume (ENS01). 

Target dwarf halos were also found in the SGIF-128 
simulation. However, because of their extremely low mass 
(corresponding to 5-6 particles in SGIF-128), it was neces- 
sary to create a second "parent" simulation for them by 
resimulating a region of the SGIF-128 volume at signifi- 
cantly higher resolution. To this end, a spherical region of 
radius 4.4 ft -1 Mpc, with mean density close to the univer- 
sal average, was selected at random within the 35.325 h~ 
Mpc box. This spherical region was then resimulated with 
roughly one hundred times more particles than in SGIF-128. 
The target dwarf halos were identified within this spherical 
volume again from an FoF(0.2) group list. A total of eigh- 
teen halos with 450-550 particles (corresponding to masses 
of 9-11 x 10 9 ft -1 Mq) were chosen. We report results on the 
four halos in this list that have been resimulated to date. 
High resolution initial conditions for these dwarf halos were 
created in an identical way to the more massive galaxy and 
cluster halos. 

Numerical parameters were chosen to ensure that all 
halos, regardless of mass, were resimulated at comparable 
mass resolution (typically over 10 6 particles within the virial 
radius at z = 0, see Tabic [SJ. 



2.6 The Analysis 

We focus our analysis on the spherically-averaged mass pro- 
file of simulated halos at z = 0. This is measured by sorting 
particles in distance from the center of each halo and ar- 
ranging them in bins of equal logarithmic width in radius. 
Density profiles, p(r), are computed simply by dividing the 
mass in each bin by its volume. The cumulative mass within 
each bin, M(r), is then used to compute the circular veloc- 
ity profile of each halo, V c (r) = W GM(r) /r, as well as the 
cumulative density profile, p{r) = 3M(r)/4nr 3 , which we 
shall use in our analysis. 

The center of each halo is determined using an iter- 
ative technique in which the center of mass of particles 
within a shrinking sphere is computed recursively until a few 
thousand particles are left (see P03 for details). In a multi- 
component system, such as a dark halo with substructure, 
this centering procedure isolates the densest region within 
the largest subcomponent. In more regular systems, the cen- 
ter so obtained approximately coincides with the centers de- 
fined by the center of mass weighted by the local density or 
by the gravitational potential of each particle. 

We note that, unlike in NFW, no attempt has been 
made to select halos at a particularly quiet stage in their 
dynamical evolution; our sample thus contains halos in equi- 
librium as well as a few with prominent substructure as a 
result of recent accretion events. 



2.7 Parameter selection criteria 

The analysis presented in P03 demonstrated that the mass 
profile of a simulated halo is numerically robust down to 
a "convergence radius", r conv , that depends primarily on 
the number of particles and time steps, as well as on the 
choice of gravitational softening in the simulation. Each of 
these choices imposes a minimum radius for convergence, al- 
though for an "optimal" choice of parameters (i.e., one that, 
for given r conv , minimizes the number of force computations 
and time steps) the most stringent criterion is that imposed 
by the number of particles within r2oo- In this optimal case, 
the minimum resolved radius is well approximated by the lo- 
cation at which the two-body relaxation time, t ro iax, equals 
the age of the universe (see Hayashi et al 2003 for further 
validation of these criteria, but see Binney 2003 for a differ- 
ent opinion). 

To be precise, we shall identify r conv with the radius 
where i ro i a x equals the circular orbital timescale at the virial 
radius, t cilc (r2oo) = 27r r 200/^200- Thus, r conv is defined by 
the following equation, 

trdax(r) _ N rjVc _ 1 _ 
£circ(r 2 oo) 8 In A r 2 oo/V2oo 

Here N = N(r) is the number of particles enclosed within r, 
and V200 = V c (r2oo). With this definition, the convergence 
radius in our best-resolved halos, outside which V c (r) con- 
verges to better than 10%, is of order ~ 0.005 r2oo- 



3 RESULTS 

3.1 Density Profiles 

The top panels of Figure show the density profiles, p(r), 
of the nineteen simulated halos in our sample. In physical 
units, the profiles split naturally into three groups: from left 
to right, "dwarf" (dotted), "galaxy" (dashed), and "cluster" 
(dot-dashed) halos, respectively. Each profile is shown from 
the virial radius, r^oo, down to the innermost converged ra- 
dius, r C onvi a convention that we shall adopt in all figures 
throughout this paper. 

The thick solid lines in the top-left panel show the NFW 
profiles (eq.0 expected for halos in each group, with param- 
eters chosen according to the prescription of Eke, Navarro 
& Steinmetz (2001). Note that these NFW curves are not 
best fits to any of the simulations, but that they still capture 
reasonably well the shape and normalization of the density 
profiles of the simulated halos. 

The top right panel of Figure is similar to the top 
left one, but the comparison is made here with the modified 
form of the NFW profile proposed by M99 (eq. 0. There 
is no published prescription specifying how to compute the 
numerical parameters of this formula for halos of given mass, 
so the three profiles shown in this panel are just "eyeball" 
fits to one halo in each group. Like the NFW profile, the 
M99 formula also appears to describe reasonably well the 
gently-curving density profiles of ACDM halos. 

Figure thus confirms a number of important trends 
that were already evident in prior simulation work. 

• ACDM halo density profiles deviate significantly from 
simple power laws, and steepen systematically from the cen- 
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ter outwards; they are shallower than isothermal near the 
center and steeper than isothermal near the virial radius. 

• There is no indication of a well defined central "core" of 
constant density; the dark matter density keeps increasing 
all the way in, down to the innermost resolved radius. 

• Simple formulae such as the NFW profile (eq.0 or the 
M99 formula (eq. appear to describe the mass profile of 
all halos reasonably well, irrespective of mass, signaling a 
"universal" profile shape. Properly scaled, a dwarf galaxy 
halo is almost indistinguishable from a galaxy cluster halo. 

We elaborate further on each of these conclusions in 
what follows. 

3.1.1 NFW vsM99 fits 

Are the density profiles of ACDM halos described better by 
the NFW formula (eq.0 or by the modification proposed by 
M99 (eq.|2l? The answer may be seen in the bottom panels 
of Figure These panels show the deviations (simulation 
minus fit) from the best fits to the density profiles of each 
halo using the NFW profile or the M99 profile. These fits 
are obtained by straightforward x 2 minimization, assigning 
equal weight to each radial bin. This is done because the 
statistical (Poisson) uncertainty in the determination of the 
density within each bin is negligible (each bin contains from 
several thousand to several hundred thousand particles) so 
the remaining uncertainties are likely to be dominated by 
systematics, such as the presence of substructure, varying 
asphericity, as well as numerical error, whose radial depen- 
dence is difficult to assess quantitatively (see P03). 

As shown in the bottom panels of Figure there is 
significant variation in the shape of the density profile from 
one halo to another. Some systems are fit better by eq. Q 
than by eq. [5] and the reverse is true in other cases. Over 
the radial range resolved by the simulations, p(r) deviates 
from the best fits by less than ~ 50%. NFW fits tend to 
underestimate the density in the inner regions of most halos; 
by up to 35% at the innermost resolved point. M99 fits, 
on the other hand, seem to do better for low mass halos, 
but tend to overestimate the density in the inner regions of 
cluster halos by up to 60%. We have explicitly checked that 
these conclusions are robust to reasonable variation in the 
binning used to construct the density profiles, as well as in 
the adopted minimization procedure. 

This level of accuracy may suffice for a number of obser- 
vational applications, with the proviso that comparisons are 
restricted to radii where numerical simulations are reliable; 
i.e., 7- conv < r < r2oo- Deviations from the best fits increase 
systematically towards the center, so it is likely that extrap- 
olations of either fitting formula to radii much smaller than 
r CO nv will incur substantial error. We discuss below (§ I3.ri|l 
possible modifications to the fitting formulae that may min- 
imize the error introduced by these extrapolations. 

3.2 Circular Velocity Profiles 

Many observations, such as disk galaxy rotation curves or 
strong gravitational lensing, are better probes of the cumula- 
tive mass distribution than of the differential density profile 
shown in Figure Q Since cumulative profiles are subject to 
different uncertainties than differential ones, it is important 



to verify that our conclusions regarding the suitability of 
the NFW or M99 fitting formulae are also applicable to the 
cumulative mass distribution of ACDM halos. 

The radial dependence of the spherically-averaged cir- 
cular velocity profile of all halos in our series is shown in 
Figure|5| As in Figure^ the thick solid curves in the top left 
(right) panel are meant to illustrate a typical NFW (M99) 
profile corresponding to dwarf, galaxy, and cluster halos, re- 
spectively. The bottom left and right panels show deviations 
from the best fit to each halo using the NFW or M99 profile, 
respectively. Both profiles reproduce the cumulative mass 
profile of the simulated halos reasonably well. The largest 
deviations seen are for the M99 fits, but they do not exceed 
25% over the radial range resolved in the simulations. NFW 
fits fare better, with deviations that do not exceed f0%. 

As with the density profiles, the deviations between sim- 
ulation and fits, although small, increase toward the center, 
suggesting that caution should be exercised when extrapo- 
lating these fitting formulae beyond the spatial region where 
they have been validated. This is important because obser- 
vational data, such as disk galaxy rotation curves, often ex- 
tend to regions inside the minimum convergence radius in 
these simulations. 



3.3 Radial dependence of logarithmic slopes 

We have noted in the previous subsections that systematic 
deviations are noticeable in both NFW and M99 fits to 
the mass profiles of simulated ACDM halos. NFW fits tend 
to underestimate the dark matter density near the center, 
whilst M99 fits tend to overestimate the circular velocity in 
the inner regions. The reason for this is that neither fitting 
formula fully captures the radial dependence of the density 
profile. We explore this in Figure |3] which shows the loga- 
rithmic slope, din p/dlnr = —f3(r), of all simulated halos, 
as a function of radius. Although there is substantial scatter 
from halo to halo, a number of trends are robustly defined. 

The first trend to note is that halo density profiles be- 
come shallower inward down to the innermost resolved ra- 
dius, r conv (the smallest radius plotted in Figure [3J. We see 
no indication for convergence to a well defined asymptotic 
value of the inner slope in our simulated halos, neither to 
the fc> = (3{r = 0) = 1 expected for the NFW profile (solid 
curves in Figure^ nor to the /3o = 1.5 expected in the case 
of M99 (dotted curves in same figure). 

The second trend is that the radial dependence of the 
logarithmic slope deviates from what is expected from either 
the NFW or the M99 fitting formulae. Near r conv , the slopes 
are significantly shallower than /3o = 1-5 (and thus in dis- 
agreement with the M99 formula) but they are also signifi- 
cantly steeper than expected from NFW fits. In quantitative 
terms, let us consider the slope well inside the characteristic 
radius, r_2 (where the slope takes the "isothermal" value 2 of 
P{r-2) = 2). For cluster halos, for example, at r — 0.1 r_2 

2 The characteristic radius, r_2, as well as the density at that 
radius, p_2 = p(r_2, can be measured directly from the sim- 
ulations, without reference to or need for any particular fitting 
formula. For the NFW profile, r_2 is equivalent to the scale ra- 
dius r s (see eq. 0. The density at r_2 is related to the NFW 
characteristic density, p s , by p_2 = p(r_2) = Ps/4. 
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(~ 50 h~ x kpc) the average slope is approximately —1.3, 
whereas the NFW formula predicts ~ —1.18 and M99 pre- 
dicts ~ —1.5. This is in agreement with the latest results of 
Fukushige, Kawai & Makino (2003) , who also report profiles 
shallower than r~ 1,J at the innermost converged radius of 
their simulations. A best-fit slope of r -1,3 was also reported 
by Moore et al (2001) for a dwarf galaxy halo (of mass simi- 
lar to the Draco dwarf spheroidal) , although that simulation 
was stopped at z — 4, and might therefore not be directly 
comparable to the results we present here. 

This discrepancy in the radial dependence of the loga- 
rithmic slope between simulations and fitting formulae is at 
the root of the different interpretations of the structure of 
the central density cusp proposed in the literature. For ex- 
ample, because profiles become shallower inward more grad- 
ually than in the NFW formula, modifications with more 
steeply divergent cusps (such as eq. tend to fit density 
profiles (but not circular velocity profiles) better in the re- 
gion interior to r_2. This is not, however, a sure indication 
of a steeper cusp. Indeed, any modification to the NFW pro- 
file that results in a more gradual change in the slope inside 
r_2 will lead to improved fits, regardless of the value of the 
asymptotic central slope. We show this explicitly below in 

§E2D 

3.4 Maximum asymptotic slope 

Conclusive proof that the central density cannot diverge as 
steeply as j3o = 1.5 is provided by the total mass inside 
the innermost resolved radius, r conv . This is because, at any 
radius r, the mean density, p(r), together with the local 
density, p(r), provide a robust upper limit to the asymptotic 
inner slope. This is given by /3 ma *(r) = 3(1— p(r)/p(r)) > f3 , 
under the plausible assumption that j3 is monotonic with 
radius. 

Figure 01 shows /3 max as a function of radius; clearly, ex- 
cept for possibly one dwarf system, no simulated halo has 
enough dark mass within r conv to support cusps as steep as 
r -1 ' 5 . The NFW asymptotic slope, corresponding to 0o = 1, 
is still consistent with the simulation data, but the actual 
central value of the slope may very well be shallower. We 
emphasize again that there is no indication for convergence 
to a well defined value of /3q: density profiles become shal- 
lower inward down to the smallest resolved radius in the 
simulations. 



3.5 A "universal" density profile 

Figure |3 shows also that there is a well-defined trend with 
mass in the slope of the density profile measured at r conv ~ 
0.005 to 0.01r2oo (the innermost point plotted for each pro- 
file): /3(r conv ) ~ 1.1 for clusters, ~ 1.2 for galaxies, and 
~ 1.35 for dwarfs. A similar trend was noted by Jing & Suto 
(2000), who used it to argue against a "universal" density 
profile shape. However, as discussed by Klypin et al (2001), 
this is just a reflection of the trend between the concentra- 
tion of a halo and its mass. It does not indicate any departure 
from similarity in the profile shape. Indeed, one does not ex- 
pect the profiles of halos of widely different mass, such as 
those in our series, to have similar slopes at a constant frac- 
tion of the virial radius. Rather, if the density profiles are 



truly self-similar, slopes ought to coincide at fixed fractions 
of a mass-independent radial scale, such as r_2. 

Figure shows the striking similarity between the 
structure of halos of different mass when all density pro- 
files are scaled to r_2 and p_2 = p(r~2)- The density profile 
of a dwarf galaxy halo then differs very little from that of 
a galaxy cluster 10 5 times more massive. This demonstrates 
that spherically-averaged density profiles are approximately 
"universal" in shape; rarely do individual density profiles 
deviate from the scaled average by more than ~ 50%. 

In the scaled units of Figure the NFW and M99 pro- 
files are fixed, and are shown as solid and dotted curves, 
respectively. With this scaling, differences between density 
profiles are more evident than when best fits are compared, 
since the latter — by definition — minimize the deviations. 
In Figure f° r example, it is easier to recognize the "ex- 
cess" of dark mass inside r_2 relative to the NFW profile 
that authors such as M99 and Fukushige & Makino (1997, 
2001, 2003) have (erroneously) interpreted as implying a 
steeply divergent density cusp. 

The similarity in mass profile shapes is also clear in 
Figure [HJj, which shows the circular velocity curves of all 
halos in our series, scaled to the maximum, V ma x and to 
the radius where it is reached, r max . NFW and M99 are 
again fixed curves in these scaled units. This comparison is 
more relevant to observational interpretation, since rotation 
curve, stellar dynamical, and lensing tracers are all more 
directly related to V c (r) than to p{r). Because of the reduced 
dynamic range of the y-axis, the scatter in mass profiles from 
halo to halo is more clearly apparent in the V c profiles; the 
NFW and M99 profiles appear to approximately bracket the 
extremes in the mass profile shapes of simulated halos. We 
discuss below a simple fitting formula that, with the aid of 
an extra parameter, is able to account for the variety of mass 
profile shapes better that either the NFW or M99 formula. 

3.6 An improved fitting formula 

Although the discussion in the previous subsections has con- 
centrated on global deviations from simple fitting formu- 
lae such as NFW or M99, it is important to emphasize 
again that such deviations, although significant, are actually 
rather small. As shown in Figure^ best NFW fits reproduce 
the circular velocity profiles to an accuracy of better than 
~ 10% down to roughly 0.5% of r2oo- Although this level 
of accuracy may suffice for some observational applications, 
the fact that deviations increase inward and are maximal 
at the innermost converged point suggests the desirability 
of a new fitting formula better suited for extrapolation to 
regions beyond those probed reliably by simulations. 

An improved fitting formula ought to reproduce: (i) the 
more gradual shallowing of the density profile towards the 
center; (ii) the apparent lack of evidence for convergence 
to a well-defined central power-law; and (iii) the significant 
scatter in profile shape from halo to halo. After some exper- 
imentation, we have found that a density profile where f3(r) 
is a power-law of radius is a reasonable compromise that 
satisfies these constraints whilst retaining simplicity; 

fi a {r) = -dlnp/dlnr = 2 (r/r_ 2 ) a , (4) 

which corresponds to a density profile of the form, 
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In( Po /p_ a ) = (-2/a)[(r/r_ a )°'-l]. (5) 

This profile has finite total mass (the density cuts off expo- 
nentially at large radius) and has a logarithmic slope that 
decreases inward more gradually than the NFW or M99 pro- 
file. The thick dot-dashed curves in Figures |3] and |3]show 
that eq.|S](with a ~ 0.17) does indeed reproduce fairly well 
the radial dependence of f3(r) and j3 max {r) in simulated ha- 
los. 

Furthermore, adjusting the parameter a allows the pro- 
file to be tailored to each individual halo, resulting in im- 
proved fits. Indeed, as shown in Figure |S| eq. reproduces 
the density profile of individual halos to better than ~ 10% 
over the reliably resolved radial range, and that there is no 
discernible radial trend in the residuals. This is a significant 
improvement over NFW or M99 fits, where the maximum 
deviations were found at the innermost resolved radius. The 
best-fit values of a (in the range 0.1 - 0.2) show no obvious 
dependence on halo mass, and are listed in Table |3] The av- 
erage q is 0.172 and the dispersion about the mean is 0.032. 

We note that the p a profile is not formally diver- 
gent, and converges to a finite density at the center, p = 
e 2/Q p_2 ~ 6 x 10 5 p-2 (for a = 0.15). It is unclear at this 
point whether such asymptotic behavior is a true property of 
ACDM halos or simply an artifact of the fitting formula that 
results from choosing pa = in eq. [I] The simulations show 
no evidence for convergence to a well-defined central value 
for the density, but even in the best-resolved cases they only 
probe regions where densities do not exceed ~ 10 2 p~2- This 
is, for a in the range 0.1 - 0.2, several orders of magnitude 
below the maximum theoretical limit in eq. 

We note as well that the convergence to /3q — is quite 
slow for the values of a favored by our fits. Indeed, for a — 
0.1, the logarithmic slope only reaches a value significantly 
shallower than the NFW asymptotic slope at radii that are 
well inside the convergence radius of our simulations; for ex- 
ample, fi a {r) only reaches 0.5 at r = 9.5 x 10 _7 r_2, corre- 
sponding to r ~ 0.01 pc for galaxy-sized halos. This implies 
that the p a profile is in practice "cuspy" for most astrophys- 
ical applications. Establishing conclusively whether ACDM 
halos actually have divergent inner density cusps is a task 
that awaits simulations with much improved resolution than 
those presented here. 

3.7 Comparison between fitting formulae 

Figure [7| compares the density and circular velocity profiles 
implied by the p a formula (eq. EJ with the NFW and M99 
profiles (left panels) , as well as with the fitting formula pro- 
posed by Stoehr et al (2002, hereafter SWTS) to describe 
the structure of substructure halos (right panels). 

The top left panel of Figure [7| shows that, despite its 
finite central density, the p a profile can approximate fairly 
well both an NFW profile (for a ~ 0.2) and an M99 profile 
(for a ~ 0.1) for over three decades in radius. The circular 
velocity profile for a — 0.2 is likewise quite similar to NFWs 
(bottom left panel of Figure , but the similarity to the 
shape of the M99 V c profile is less for all values of a. 

Interestingly, the V c profiles corresponding to p a resem- 
ble parabolae in a log-log plot, and thus may be used to 
approximate as well the mass profiles of substructure halos, 
as discussed by SWTS. This is demonstrated in the bottom 



right panel of Figure |7| where we show that the V c profiles 
corresponding to a = 0.1, 0.2, and 0.7, are very well approx- 
imated by the SWTS formula, 

log(V c /VUx) = ~a[log(r/r max )] 2 , (6) 

for a — 0.09, 0.17, and 0.45, respectively. The latter value 
(a — 0.45, or a — 0.7) corresponds to the median of the best 
SWTS fits to the mass profile of substructure halos. Note 
that this is quite different from the a ~ 0.1 - 0.2 required to 
fit isolated ACDM halos (see Table El . 

It might actually be preferable to adopt the p a profile 
rather than the SWTS formula for describing substructure 
halos, since p a (r) is monotonic with radius and extends over 
all space. This is not the case for SWTS, as shown in the 
top-right panel of Figure |7| The SWTS density profiles are 
"hollow" (i.e., the density has a minimum at the center), and 
extend out to a maximum radius, given by e 1//4a r max . This 
is because the circular velocity in the outer regions of the 
SWTS formula fall off faster than Keplerian, and therefore 
the corresponding density becomes formally negative at a 
finite radius. 

The p a profile thus appears versatile enough to repro- 
duce, with a single fitting parameter, the structure of ACDM 
halos and of their substructure. Since p a captures the inner 
slopes better than either the NFW or M99 profile, it is also 
likely to be a safer choice should extrapolation of the mass 
profile beyond the converged radius prove necessary. We end 
by emphasizing, however, that all simple fitting formulae 
have shortcomings, and that direct comparison with simula- 
tions rather than with fitting formulae should be attempted 
whenever possible. 

3.8 Scaling parameters 

The application of fitting formulae such as the one described 
above requires a procedure for calculating the characteristic 
scaling parameters for a given halo mass, once the power 
spectrum and cosmological parameters are specified. NFW 
developed a simple procedure for calculating the parameters 
corresponding to halos of a given mass. Because of the close 
relationship between the scale radius, r s , and characteristic 
density, p s , of the NFW profile and the r_2 and p_2 pa- 
rameters of eq. |K| we can use the formalism developed by 
NFW to compute the expected values of these parameters 
in a given cosmological model. 

NFW interpreted the characteristic density of a halo as 
reflecting the density of the universe at a suitably defined 
time of collapse. Their formalism assigns to each halo of 
mass M (identified at z — 0) a collapse redshift, z co u(M, /) 
defined as the epoch when half the mass of the halo was 
first contained in progenitors more massive than a certain 
fraction / of the final mass. With this definition, and once 
/ has been chosen, z co n can be computed using the Press- 
Schechter theory (e.g., Lacey & Cole 1993). The NFW model 
then assumes that the characteristic density of a halo (i.e., p s 
in eq. is proportional to the mean density of the universe 
at z C oii- 

The redshift dependence of the characteristic density 
was first probed in detail by Bullock et al (2001, here- 
after B01), who proposed a modification to NFWs model 
in which, for a given halo mass, the scale radius, r s , remains 
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approximately constant with redshift. Eke, Navarro & Stein- 
metz (2001, hereafter ENS), on the other hand, argued that 
the characteristic density of a halo is determined by the 
amplitude and shape of the power spectrum, as well as by 
the universal expansion history. Their formalism reproduces 
nicely the original results of NFW as well as the redshift 
dependence pointed out by B01, and is applicable to more 
general forms of the power spectrum, including the "trun- 
cated" power spectra expected in scenarios such as warm 
dark matter (see ENS for more details). 

We have used the ENS and B01 formalisms to predict 
the halo mass dependence of the scaling parameters, p-2 
and r_2, and we compare the results with our simulations 
in Figure [S] The ENS prediction is shown by the solid line 
whereas the dotted line shows that of B01. Both formalisms 
reproduce reasonably well the trend seen in the simulations, 
so that one can use either, in conjunction with eq. |K| (with 
a in the range 0.1-0.2), to predict the structure of a ACDM 
halo. A simple code that computes r_2 and p_2 as a function 
of mass in various cosmological models is available upon re- 
quest from the authors. Existing codes that compute NFW 
halo parameters as a function of mass and of other cosmo- 
logical parameters may also be used, noting that p_2 = Ps/4 
and that r_2 = r s . 

Finally, we note that neither formalism captures per- 
fectly the mass dependence of the characteristic density; 
small but significant deviations, as well as a sizable scatter, 
are evident in Figure |H| Dwarf galaxy halos appear to be 
less concentrated than predicted by the formalism proposed 
by B01; a similar observation applies to cluster halos when 
compared to ENS' predictions. Such shortcomings should be 
considered when deriving cosmological constraints from fits 
to observational data (see, e.g., Zentner & Bullock 2002, Mc- 
Gaugh et al 2003); and suggest again that direct comparison 
between observation and simulations is preferable to the use 
of fitting formulae. 



4 SUMMARY 

We have analyzed the mass profile of ACDM halos in a series 
of simulations of high mass, spatial, and temporal resolu- 
tion. Our series targets halos spanning five decades in mass: 
"dwarf" galaxy halos with virial circular velocities of order 
V200 ~ 30 km s _1 ; "galaxy" -sized halos with V200 ~ 200 km 
s _1 ; and "cluster" halos with V200 ~ 1200 km s _1 . Each of 
the nineteen halos in our series was simulated with com- 
parable numerical resolution: they have between 8 x 10 5 
and 4 x 10 6 million particles within the virial radius, and 
have been simulated following the "optimal" prescription 
for time-stepping and gravitational softening laid down in 
the numerical convergence study of P03. 

The high resolution of our simulations allows us to 
probe the inner properties of the mass profiles of ACDM 
halos, down to ~ 0.5% of r2oo in our best resolved runs. 
These results have important implications for the structure 
of the inner cusp in the density profile and resolve some of 
the disagreements arising from earlier simulation work. Our 
main conclusions may be summarized as follows. 

• ACDM halo density profiles are "universal" in shape: 
i.e., a simple fitting formula reproduces the structure of all 
simulated halos, regardless of mass. Both the NFW profile 



and the profile proposed by M99 describe the density and 
circular velocity profiles of simulated halos reasonably well. 
Best NFW fits to the circular velocity profiles deviate by 
less than 10% over the region which is well resolved numer- 
ically. Best M99 fits reproduce circular velocity profiles to 
better than 25% over the same region. It should be noted, 
however, that the deviations increase inwards and are typ- 
ically maximal at the innermost resolved radius, a result 
that warns against extrapolating to smaller radii with these 
fitting formulae. 

• ACDM halos appear to be "cuspy": i.e., the dark mat- 
ter density increases monotonically towards the center with 
no evidence for a well-defined "core" of constant density. We 
find no evidence, however, for a central asymptotic power- 
law bin the density profiles. These become progressively 
shallow inwards and are significantly shallower than isother- 
mal at the innermost resolved radius, r conv . At r ~ 0.01 V2oo, 
the average slope of "cluster", "galaxy" and "dwarf" halos 
halos is /3(r conv ) ~ 1.1, ~ 1.2, and ~ 1.35, respectively. This 
is steeper than predicted by the NFW profile but shallower 
than the asymptotic slope of the M99 profile. 

• The density and enclosed mass at r conv may be used to 
derive an upper limit on any asymptotic value of the inner 
slope. Cusps as steep as f3o = 1.5 are confidently ruled out in 
essentially all cases; the asymptotic slope of the NFW profile 
(/So = 1) is still consistent with our data. The radial depen- 
dence of f3(r) differs from that of the NFW profile, however, 
decreasing more slowly with decreasing radius than is pre- 
dicted. For some scalings of the NFW fitting formula to the 
numerical data, this shape difference appears as a dark mat- 
ter "excess" near the center which has (erroneously) been 
interpreted indicating a steeply divergent density cusp. 

• A simple formula where /3(r) is a power law of radius 
reproduces the gradual radial variation of the logarithmic 
slope and its apparent failure to converge to any specific 
asymptotic value (eq. 0. This formula leads to much im- 
proved fits to the density profiles of simulated halos, and 
may prove a safer choice when comparison with observation 
demands extrapolation below the innermost converged radii 
of the simulations. 

Our study demonstrates that, although simple fitting 
formulae such as NFW are quite accurate in describing the 
global structure of ACDM halos, one should be aware of 
the limitations of these formulae when interpreting obser- 
vational constraints. Extrapolation beyond the radial range 
where these formulae have been validated is likely to produce 
substantial errors. Proper account of the substantial scatter 
in halo properties at a given halo mass also appears neces- 
sary when assessing the consistency of observations with a 
particular cosmological model. Direct comparison between 
observations and simulations (rather than with fitting for- 
mulae) is clearly preferable whenever possible. Given the 
computational challenge involved in providing consistent, ro- 
bust, and reproducible theoretical predictions for the inner 
structure of CDM halos it is likely that observational con- 
straints will exercise to the limit our hardware and software 
capabilities for some time to come. 
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Table 1. Parameters of the parent cosmological simulations 

Label ^box 

[ft -1 Mpcl 



[h 



ENS01 32.5 
SGIF-128 35.325 
ACDM-512 479.0 



49.0 1.36 x 10 9 
49.0 1.75 x 10 9 
36.0 6.82 x 10 10 



10 
10 

30 



CODE 



AP3M 
GADGET 
GADGET 



Table 2. Main parameters of rcsimulatcd halos 



Label 


Zi 


e 

[h- 1 kpc] 


^200 


M 2 oo 
[h- 1 M ] 


T200 
[h- 1 kpc] 


V200 
[km s- 1 ] 


^conv 

[h- 1 kpc] 


CODE 


Dl 


74 


0.0625 


784980 


7.81 x 10 9 


32.3 


32.3 


0.34 


GADGET 


D2 


49 


0.0625 


778097 


9.21 x 10 9 


34.1 


34.1 


0.37 


GADGET 


D3 


49 


0.0625 


946421 


7.86 x 10 9 


32.3 


32.3 


0.33 


GADGET 


D4 


49 


0.0625 


1002098 


9.72 x 10 9 


34.7 


34.7 


0.32 


GADGET 


Gl 


49 


0.15625 


3447447 


2.29 x 10 12 


214.4 


214.4 


1.42 


GADGET 


G2 


49 


0.5 


4523986 


2.93 x 10 12 


232.6 


232.6 


1.25 


PKDGRAV 


G3 


49 


0.45 


2661091 


2.24 x 10 12 


212.7 


212.7 


1.65 


PKDGRAV 


G4 


49 


0.3 


3456221 


1.03 x 10 12 


164.0 


164.0 


1.01 


PKDGRAV 


G5 


49 


0.35 


3913956 


1.05 x 10 12 


165.0 


165.0 


1.02 


PKDGRAV 


G6 


49 


0.35 


3739913 


9.99 x 10 11 


162.5 


162.5 


1.03 


PKDGRAV 


G7 


49 


0.35 


3585676 


9.58 x 10 11 


160.3 


160.3 


1.02 


PKDGRAV 


CI 


36 


5.0 


1565576 


7.88 x 10 14 


1502.1 


1502.1 


16.8 


GADGET 


C2 


36 


5.0 


1461017 


7.36 x 10 14 


1468.1 


1468.1 


16.9 


GADGET 


C3 


36 


5.0 


1011918 


5.12 x 10 14 


1300.6 


1300.6 


16.1 


GADGET 


C4 


36 


5.0 


1050402 


5.31 x 10 14 


1316.7 


1316.7 


15.9 


GADGET 


C5 


36 


5.0 


1199299 


6.05 x 10 14 


1375.5 


1375.5 


16.2 


GADGET 


C6 


36 


5.0 


1626161 


8.19 x 10 14 


1521.1 


1521.1 


15.5 


GADGET 


C7 


36 


5.0 


887837 


4.50 x 10 14 


1245.8 


1245.8 


16.4 


GADGET 


C8 


36 


5.0 


1172850 


5.92 x 10 14 


1365.4 


1365.4 


16.8 


GADGET 



Table 3. Fit and structural parameters of resimulatcd halos 

Label r_2 p-2 r max V ma x 



[h- 1 kpc] 



P-2 
[Pcrit] 



'"max 

[h- 1 kpc] 



[ku 



[h- 1 kpc] 



[Pcrit] 



[h- 1 kpc] 



PM 
[Pcrit] 



Dl 


3.23 


1.12e4 


6.07 


39.1 


2.59 


7.03e4 


5.38 


7.58e3 


0.164 


D2 


3.04 


1.58e4 


6.64 


44.2 


2.43 


9.61e4 


2.27 


1.17e5 


0.211 


D3 


2.57 


1.58e4 


6.35 


36.9 


2.94 


5.01e4 


4.05 


2.51e4 


0.122 


D4 


2.57 


2.24e4 


4.27 


45.7 


2.06 


1.49e5 


2.18 


1.36e5 


0.166 


Gl 


18.5 


6.76e3 


23.7 


1.95e2 


23.2 


4.06e4 


19.4 


8.43e4 


0.142 


G2 


28.0 


2.40e3 


68.5 


1.78e2 


16.8 


1.13e5 


19.4 


8.43e4 


0.191 


G3 


20.2 


6.31e3 


43.4 


1.96e2 


28.0 


1.52e4 


47.3 


8.24e3 


0.142 


G4 


29.6 


4.37e3 


63.4 


2.49e2 


12.3 


6.78e4 


16.8 


3.44e4 


0.177 


G5 


20.7 


1.58e4 


67.7 


2.91e2 


13.8 


5.20e4 


15.3 


4.23e4 


0.184 


G6 


39.6 


2.00e3 


96.4 


2.26e2 


15.3 


3.79e4 


20.7 


2.03e4 


0.171 


G7 


16.4 


1.26e4 


29.9 


1.94e2 


13.4 


6.22e4 


14.9 


5.15e4 


0.138 


CI 


5.84e2 


4.68e2 


1.03e3 


1.48e3 


440 


3.36e3 


661 


1.58e3 


0.133 


C2 


3.95e2 


1.15e3 


9.99e2 


1.51e3 


362 


5.17e3 


396 


4.46e3 


0.215 


C3 


3.27e2 


1.12e3 


6.15e2 


1.38e3 


249 


9.07e3 


278 


7.44e3 


0.188 


C4 


4.16e2 


7.94e2 


6.57e2 


1.38e3 


315 


5.47e3 


339 


4.91e3 


0.161 


C5 


2.87e2 


1.91e3 


6.48e2 


1.42e3 


271 


8.45e3 


326 


5.88e3 


0.215 


C6 


3.82e2 


1.32e3 


6.94e2 


1.64e3 


297 


8.70e3 


302 


8.75e3 


0.203 


C7 


5.69e2 


3.55e2 


1.25e3 


1.25e3 


283 


3.92e3 


475 


2.11e3 


0.129 


C8 


3.68e2 


1.00e3 


9.35e2 


1.44e3 


361 


4.41e3 


345 


5.04e3 


0.219 
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Figure 1. Spherically-averaged density profiles of all our simulated halos. Densities are computed in radial bins of equal logarithmic 
width and are shown from the innermost converged radius (r conv ) out to about the virial radius of each halo (r2uo)- Our simulations 
target halos in three distinct mass groups: "dwarf", "galaxy", and "cluster" halos. These groups span more than five decades in mass. 
Thick solid lines in the top panels illustrate the expected halo profile for each mass range according to the fitting formula proposed by 
NFW (top-left) or M99 (top-right). Bottom panels indicate the deviation from the best fit achieved for each individual halo (simulation 
minus fit) with the NFW profile (eq. or with its modified form, as proposed by M99 (eg. 121. 
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Figure 2. Spherically-averaged circular velocity (V c (r) = y/ GM(r)/r) profiles of all our simulated halos. As in Figure circular 
velocities are computed in radial bins of equal logarithmic width and are shown from the innermost converged radius (r con v) out to 
about the virial radius (r20o) of each halo. Our simulations target halos in three distinct mass groups: "dwarf", "galaxy", and "cluster" 
halos, spanning more than a factor of ~ 50 in velocity. Thick solid lines in the top panels illustrate the expected profile for each mass 
range according to the fitting formula proposed by NFW (top-left) or M99 (top-right). Bottom panels indicate the deviation from the 
best fit achieved for each individual halo (simulation minus fit) with the NFW profile (eq. or with its modified form, as proposed by 
M99 (eq.0. 
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Figure 3. Logarithmic slope of the density profile of all halos in our sample, plotted versus radius. Thick solid and dotted curves 
illustrate the radial dependence of the slope expected from the NFW profile (eq. and the modification proposed by M99 (eq. 151 . 
respectively. Note that although both fitting formulae have well-defined asymptotic inner slopes (—1 and —1.5, respectively) there is no 
sign of convergence to a well-defined value of the central slope in the simulated halos. At the innermost converged radius, the simulated 
halo profiles are shallower than —1.5, in disagreement with the Moore et al profile. Also, inside the radius at which the slope equals —2, 
r_2, the profiles appear to get shallower more gradually than in the NFW formula. A power-law radial dependence of the slope seems 
to fit the results of our simulations better; the dot-dashed lines indicate the predictions of the p a profile introduced in eqs, HI and 151 for 
a = 0.17. Best fits to individual halos yield a in the range 0.1-0.2 (see Table 131 . 
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Figure 4. Maximum asymptotic inner slope compatible with the mean density interior to radius r, p(r), and with the local density at 
that radius, p(r). This provides a robust limit to the central slope, fSo < /3 ma x(r) = —3(1 — p(r)/p(r)), under the plausible assumption 
that /3 is monotonic with radius. Note that there is not enough mass within the innermost converged radius in our simulations to support 
density cusps as steep as r~ 1,5 . The asymptotic slope of the NFW profile, f3o = 1, is still compatible with the simulated halos, although 
there is no convincing evidence for convergence to a well defined power-law behavior in any of our simulated halos. The thick dot-dashed 
curves illustrate the expected radial dependence of /3 max for the p a profile introduced in § 13.61 for a = 0.17. 
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Figure 5. (a-left) Density profiles of all halos in our series, scaled to the radius, r_2, where the local logarithmic slope of the density 
profile takes the isothermal value of /3 = — d log p/d log r = 2. Densities are scaled to p_2 = p(r_2). This figure shows that, with proper 
scaling, there is little difference in the shape of the density profile of halos of different mass, confirming the "universal" nature of the 
mass profile of ACDM halos. The NFW profile (eq. is a fixed curve in these scaled units, and is shown with a thick solid line. The M99 
formula (eq. [3J is shown with a dashed line, (b-right) Circular velocity profiles all halos in our series, scaled to the maximum velocity, 
Vmax, and to the radius at which it is reached, r max . Note the significant scatter from halo to halo, and also that the NFW and M99 
profiles appear to bracket the extremes of the mass profile shapes of halos in our simulation series. 
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Figure 6. As Figure^ but for the p a fitting formula presented in eq.|2] Thick solid lines in the top panels illustrate the expected halo 
profile for each mass range according to the prescription proposed by NFW. Bottom panels indicate the deviation from the best p a fit 
achieved for each individual halo, taking a as a free parameter. Note the improvement in the fits compared with those achieved with the 
NFW or M99 profile and shown in Figure 
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Figure 7. Comparison between the density (top) and circular velocity (bottom) profiles corresponding to four different fitting formulae: 
NFW (solid curves, eq. □ , M99 (dotted curves, eq. 0, SWTS (dashed curves, eq. and p a (dot-dashed curves, eq. Circular velocity 
profiles are scaled to the maximum, V ma x, and to the radius where that is reached, r max . Density profiles are scaled as in Figure IB1 
Note that, despite having a finite central density, the p a formula matches, for about 3 decades in radius, the NFW profile (for a = 0.2) 
or the M99 profile (for a = 0.1, see top left panel). It also matches closely the SWTS "parabolic" circular velocity profiles intended to 
reproduce substructure halos (see bottom right panel); the V c profile with a = 0.7 is very similar to the SWTS profile with a = 0.45, the 
median value of the fits to substructure halos reported by SWTS. See text for further discussion. 
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Figure 8. The radius, r_2, where the logarithmic slope of the density profile takes the "isothermal" value, f)(r-2) = 2, plotted versus 
the local density at that radius, p_2 = p(r_2), for all simulated halos in our series. This figure illustrates the mass dependence of 
the central concentration of dark matter halos: low mass halos are systematically denser than their more massive counterparts. Solid 
and dotted lines indicate the scale radius-characteristic density correlation predicted by the formalisms presented by Eke, Navarro & 
Steinmetz (2001) and Bullock et al (2001). These parameters may be used, in conjunction with eq. |5] to predict the mass profile of 
ACDM halos. 



